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We calculate the mean and variance of net-baryon number and net-electric charge distributions 
from Quantum Chromodynamics (QCD) using a next-to-leading order Taylor expansion in terms 
of temperature and chemical potentials. We compare these expansions with experimental data 
from STAR and PHENIX, determine the freeze-out temperature in the limit of vanishing baryon 
chemical potential, and, for the first time, constrain the curvature of the freeze-out line through 
a direct comparison between experimental data on net-charge fluctuations and a QCD calculation. 

We obtain a bound on the curvature coefficient, < 0.011, that is compatible with lattice QCD 
results on the curvature of the QCD transition line. 

PACS numbers: ll.15.Ha, 12.38.Gc, 12.38.Mh, 24.60.-k 


INTRODUCTION 

Heavy ion collisions at varying beam energies are per¬ 
formed at the Relativistic Heavy Ion Collider (RHIC) 
with the goal to probe properties of strong-interaction 
matter at different temperatures and baryon chemical 
potentials. This Beam Energy Scan program aims for ex¬ 
ploring physics in the vicinity of the pseudo-critical line 
of the transition from hadronic matter to the quark-gluon 
plasma. The hope is to find evidence for the existence 
of a critical point [l[ in the phase diagram of strong- 
interaction matter that marks the end of a line of first 
order phase transitions. 

With decreasing beam energy, s/sffff, more baryons 
from the incident nuclei are stopped leading to an in¬ 
crease in the baryon number density and thus also to 
an increase of the baryon chemical potential in a cen¬ 
tral rapidity window. The dense matter created in these 
collisions is expected to reach local thermal equilibrium 
quickly. It subsequently expands and cools down. Even¬ 
tually hadrons start to form again. This hadronization is 
characterized by a temperature (T) and baryon chemical 
potential (/is) that varies with Shortly thereafter 

inelastic interactions among the hadrons cease to occur 
and different particle species ” freeze-out”. This so-called 
chemical freeze-out is characterized by a set of freeze- 
out parameters, (Ty, yi^). With varying y/sffff they map 
out a line in the T-^lb plane, called the freeze-out line, 
Tf{iiB)- An important question is how close the freeze- 
out line is to the pseudo-critical line, Tc^Rb) ii , char¬ 
acterizing the crossover transition of QCD. The proxim¬ 
ity of both lines in the phase diagram is a pre-requisite for 


being able to explore critical behavior in the vicinity of a 
possibly existing critical point, by analyzing observables 
that can probe physics on the freeze-out line. 

Unlike the crossover line, Tc{ub), for which the tran¬ 
sition temperature at rb = 0 IJ, l5| and the leading or¬ 
der correction at small r^b 0 0| have been determined 
in lattice QCD calculations, the parametrization of the 
freeze-out line, Tf{RB), is not given in terms of funda¬ 
mental parameters of QCD. The line characterizes the ex¬ 
panding medium formed in heavy ion collisions. Freeze- 
out parameters have been determined in a wide range 
of y/sffff by comparing experimental data on particle 
yields with statistical hadronization models (HRG mod¬ 
els). Parametrizations of the freeze-out line Tf^RB) have 
been extracted from such analyses [1, [3|- Physically mo¬ 
tivated ansatze, such as identifying Tf{RB) with a line of 
constant energy per particle , naturally lead to the ex¬ 
pectation that Tj^rb) is a function of r% that decreases 
with increasing rb- In fact, a simple ansatz [^, 

Tulips) = Tffi ^I — K 2 MB ~ i^{pb^ j (1) 

with Rb = Rs/Tffi, provides a good parametrization in 
the entire energy range probed in heavy ion collisions. 
Using this ansatz for the comparison of HRG model cal¬ 
culations with experimental data obtained for a wide 
range of energies ^/s^ gave = 0.023(3) 0 for the 
curvature coefficient. However, a recently performed re¬ 
fined hadronization model analysis suggests a weaker en¬ 
ergy dependence of the freeze-out line 0. Moreover, 
the attempt to capture more accurately the behavior 
of freeze-out temperatures at large i.e. small 

RB, lead to a phenomenological parametrization 0 that 
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does not even have a power-like dependence on he. 
Tf{iJ.B) — Tf^o ^ exp(—a/zTs) which favors ~ 0. 

In this paper we improve over the current situation 
by outlining a procedure to determine the 0[y?g) coef¬ 
ficient of the freeze-out line from measurements of the 
mean and variance of net-electric charge and net-proton 
number distributions and, for the first time, illustrate 
how this can be done by comparing experimental data 
obtained at large directly with a QCD calcula¬ 

tion. This puts the determination of the curvatures of 
the freeze-out line on a par with that of the QCD tran¬ 
sition line. 


RATIO ON CHARGE FLUCTUATIONS ON THE 
FREEZE-OUT LINE 


The mean, Mx = Xi (T, m) i s-^d variance, = 
X^(T, zi), of net-electric charge {X = Q) and net-baryon 
number {X = B) distributions are obtained as functions 
of T and ^ = (pB, ^^s) by taking derivatives of the 
QCD pressure with respect to pLx = l^x/T 

Xu{T,p) = ^^ , X = B, Q, S . (2) 

OPx 

The ratios of mean and variance. 





X^(T,P) 

X^iT,p) ’ 


(3) 


can be analyzed in heavy ion experiments. Implement¬ 
ing the constraints Ms = 0 and Mq/Mb = r-, which 
are appropriate for the initial conditions met in such col¬ 
lisions, the ratios and i ?^2 become functions of T 
and pb only and = ra’^jag. In lead¬ 

ing order (LO) Taylor expansion is independent of 
ji b , while the ratios and depend linearly on ps, 
R ^2 = They therefore provide mutu¬ 

ally independent information that can be used to extract 
Tf{pB) up to e>(zi|) 0. 

In order to determine the freeze-out temperature at 
Pb = 0 and the curvature of the freeze-out line, we con¬ 
sider a next-to-leading order Taylor expansion of in 
terms of T and pB around the point {Tf^^pB = 0) with 
PS and pQ being implicit functions of T and pB- Using 
Eq. [T]as a parametrization of the freeze-out line, we find 
in next-to-leading order (NLO), 






dE ?®’0 


dT 


Al • (4) 


Tf,o, 


The LO expansion coefficient is easily related to the 
quadratic fluctuations of net-electric charge and net- 
baryon number, E^^’° = rxi{T)/x 2 {T) at zero pB- 
The NLO expansion coefficient E^^’^ depends on fourth 
order cumulants, which also can be calculated in a lattice 


QCD calculation at vanishing pB- An explicit expression 
for E^®’^ will be given elsewhere [l0 . 

In order to facilitate a comparison with the experimen¬ 
tal data it is of advantage to eliminate p b from Eq. |4] in 
favor of observables that are accessible to experiments 
and QCD calculations. For a consistent treatment of the 
NLO result, Eq. |4l it suffices to use the LO relation be¬ 
tween pb and the ratio i?f 2 , i.e. pB = Bi 2 /Ri 2 ^- Tbe 
LO expansion coefficient Ri 2 ^ has been evaluated before 
and continuum extrapolated results for mi = ms/20 ob¬ 
tained with the HISQ action have been shown in [l^ . 

After replacing pB favor of the NLO Taylor 
expansion of E^-®, introduced in Eq. |4l, becomes 


E?^ = ai2 (l + ci2 (O') + O ((i?f2)") , (5) 

where ai 2 (r) = The coefficient of the quadratic 

correction, Ci 2 , depends on the parametrization of the 
freeze-out line and needs to be determined at T/q, 

C12(T/,o, 4) = Ci2{Tfp) - 4A’i2(T/,o) (6) 


with 


-12 


(r) = 


R 


Di2(T) = 


B,1 
12 , 
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yQB,2 

,0 
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B,l 
12 , 
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.dlnE^^’O 
dT 


(7) 


In Fig. [T] we show results for and 0^2 

obtained in lattice QCD calculations with the (2-I-I)- 
flavor HISQ action ll| using a physical value of the 
strange quark mass and two sets of light quark masses, 
mi = ms/20, msiXJ, which in the continuum limit corre¬ 
spond to pion mass values, ~ 160 MeV and 140 MeV, 
respectively. Here we used r ~ 0.4 which is appropriate 
for describing colliding gold or lead nuclei. Further de¬ 
tails on the simulation parameters can be found in [l0. 
We note that quark mass effects are small for the observ¬ 
ables under consideration. 

The continuum extrapolated results for E^®’° and 
i?f 2 ^ in the left and middle panels of Fig. [T] were ob¬ 
tained by performing cubic spline fits to all mi = mg/20 
data, with 1/N^ dependence for the spline coefficients 
and for the varying locations of three knots. These fits 
were performed over many bootstrap samples drawn from 
the Gaussian errors of data points and have been con¬ 
strained to agree with HRG at T = 130 MeV within 
10%. The hnal continuum results were obtained from 
mean values and errors of these bootstrapped fit results 
weighted by the quality of the fits given by the Akaike 
information criteria. The continuum extrapolations are 
consistent with our earlier results . However, statis¬ 
tical errors are reduced considerably. 

The parameter D 12 is obtained from the continuum 
extrapolated results for R^^ 


and E^^’O. It is shown 


as 
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FIG. 1. Temperature dependence of the LO expansion coefficients of the ratio (left) and R^2 (middle) and the NLO 
coefficient C 12 (right) introduced in Eq. [7] The NLO coefficient D 12 (Eq. [7| is shown as an insertion in the middle panel. 
Shown are data from calculations with the HISQ action on lattices of size x Nr, with No- = ANr- Bands show continnum 
extrapolations for the LO observables and continuum estimates for the NLO observable C 12 (see text). 


an insertion in the middle panel of Fig. [TJ The coeffi¬ 
cient C 12 receives contributions from fourth order cumu- 
lants and thus is more difficult to extract. At present we 
have calculated it for two lattice spacings and therefore 
only provide an estimate for its continuum limit. This 
is shown in the right panel of Fig. [TJ For the continuum 
estimate of 0^2 we followed an identical procedure of the 
cubic spline fits outlined above but with only two knots. 

A determination of and i?f 2 at large suffices 

to fix the freeze-out temperature Ty,o and the quadratic 
correction C 12 . Combining this with lattice QCD results 
on C 12 and D 12 allows to extract the curvature, of 
the freeze-out line. However, in heavy ion collisions 
and Rf 2 are not directly accessible. Net-proton rather 
than net baryon numbers, i.e. ratios like R 12 and E^^ 
rather than and E^^, are measured. The STAR 
Collaboration obtained Rf 2 hr the transverse momentum 
interval 0.4 GeV < Pt < with = 0.8 GeV 
(STAR0.8) 1^. The pt-range has recently been extended 
and preliminary results up to = 2.0 GeV (STAR2.0) 
[l5| have been presented. The ratio of net-electric charge 
fluctuations, R ?^, has been measured by STAR 0 and 
PHENIX 171 in the interval p™™ < Pt < 2.0 GeV with 
pmzn ^ Q 2 GeV and = 0.3 GeV[2^, respectively. 

Fig. [2] shows results for E^^ versus . In the case 

of the PHENIX data for we used the STAR2.0 data 
set to construct E^^. 

Before entering into details of the analysis of these data 
using lattice QCD calculations for i?f 2 s-nd R ?2 need 
to discuss systematic effects that arise in equilibrium 
thermodynamics because net-proton rather than net 
baryon numbers are measured in experiments and be¬ 
cause only data from a limited region in momentum space 
is available. The influence of low and high momentum 
cuts on charge fluctuations has been analyzed in HRG 


models [1^ |19| . The most important effects arise from a 


non-zero p™*" which most drastically influences the pion 


contributions to net-electric charge fluctuations. Imple¬ 
menting the pt-cuts STAR has used for protons and other 
charged particles [l^ in a HRG model calculation sug¬ 
gests that i ?^2 is overestimated by about 5% while the 
larger = 0.3 GeV used by PHENIX 171 amounts 
to an increase of 20% [l^. This explains about 40% of 
the difference in i ?^2 seen by STAR and PHENIX. It is 
conceivable that the remainder arises from the small az¬ 
imuthal coverage in the PHENIX experiment which fur¬ 
ther reduces the acceptance of charged particles [i^. In 
order to get better controll over these effects a more de¬ 
tailed experimental study of systematics of arising from 
non-zero pt cuts and cuts in the rapidity coverage will be 
needed. 


With proton number not being a conserved quantity 
[^ . the ratio clearly has no meaning in the high 
temperature plasma phase of QCD. At chemical freeze- 
out, however, when inelastic interactions are no longer 
of relevance, the net proton number may be considered 
to be a well defined concept. Still R 12 and i?f 2 will 
differ in equilibrium thermodynamics. The difference 
can be estimated in a HRG model calculation where 
= tanh(/},B-|-/I q), independent of the value ps, while 
i?f 2 explicitly depends on ps (Eq.|3|). For our NLO anal¬ 
ysis of cumulant ratios it suffices to determine the LO 
relation between and Using Taylor expansions 
for both, this obviously yields R 12 /R 12 = which is 

shown in the middle panel of Fig. [I] 

The above discussion suggests that in the STAR mea¬ 
surements the systematic errors of i ?^2 and i?f 2 are of 
similar magnitude and tend to cancel to a large extent 
in the ratio E^-®; i.e. the ratio E^-^ indeed seems to be 
a good proxy for E^-® while for the PHENIX set-up it 
overestimates E^® by at least 10%. 

We have fitted the three data sets corresponding to the 
published STAR data (STAR0.8), the preliminary STAR 
data (STAR2.0) and the PHENIX data on net-charge 
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FIG. 2. The ratio of ratios of mean and variance for net- 
electric charge and net-proton number fluctuations measured 
by the STAR and PHENIX Collaborations. For the STAR 
electric charge data we show results obtained by nor¬ 
malizing with net-proton results from the published STAR0.8 
data set [3 (triangles) as well as the preliminary STAR2.0 
data set (circles). The experimental data is compared to 
QCD predictions for two values of Kj text). The electric 
charge results obtained by PHENIX [l^ have been normal¬ 
ized by using the STAR2.0 data set on net-proton fluctuations 
(boxes). Eor orientation the upper x-axis shows ^sjviv ener¬ 
gies of the RHIC beam energy scan with labels put at the 
values for {Mp/ap)^ corresponding to the STAR2.0 data set. 
Errors on {Mp/ap)^ are not visible as they are smaller than 
the size of the symbols. 


(MQ/ag,)/(Mp/o2) 

QCD: (T,„, 4=0) 

(T,0,4=0.02) 
STAR: p™’‘=2.0 GeV . 
p,"'“‘=0.8 GeV . 
PHENIX/STAR2.0 ' 



fluctuations normalized to the STAR data on net-proton 
fluctuations (PHENIX/STAR2.0) using the ansatz given 
in Eq. [S] These data are shown in Eig. [21 For each of 
these three data sets two different fit ranges have been 
chosen, i?f 2 < 0.6 and 0 . 8 , which correspond to data 
taken at > 39 GeV and 27 GeV, respectively. This 

determines the intercept ai 2 and the curvature parameter 
ci 2 given in Table HI Differences arising from the two fit 
ranges have been added as systematic error in the error 
analysis of 012 and C 12 . From the intercept at 77^2 = 0 

the freeze-out temperature is obtained. Once Tj^ 
is fixed this way the NLO expansion coefficients c °2 and 
Di 2 are also fixed (see Table |T| and we obtain QGD pre¬ 
dictions for with as the sole free parameter. In 
Fig. [2] we show as a function of for = 0 

(colored bands). Note that any value > 0 will result in 
a weaker dependence of E^-® on i?f 2 - an illustration 
we also show the result for E^-® at Tf^o and K 2 = 0 . 02 . 
as black lines. 


CURVATURE OF THE FREEZE-OUT LINE 

We now can discuss constraints for the curvature co- 

f 

efficient K 2 resulting from the measured cumulant ratios 
E^-®. We consider the STAR data where cumulant ra¬ 
tios of net-electric charge as well as net-proton number 



STAR 0.8 

STAR2.0 

PHENIX/STAR2.0 

ai2 

0.079(3) 

0.087(2) 

0.110(9) 

C12 

0.858(101) 

0.329(74) 

0.559(352) 

T /,0 [MeV] 

145(2) 

147(2) 

155(4) 

cUTf.o) 

0.343(31) 

0.326(32) 

0.265(52) 

Di2(r/,o) 

7.04(44) 

6.62(36) 

5.27(78) 

4 

-0.073(16) 

-0.001(12) 

-0.056(67) 


TABLE I. Parameters of a quadratic fit to the STAR data on 
the ratio and the combination of PHENIX data on 
and STAR data on Ri 2 - The third to fifth row give the freeze- 
out temperature T/,o, C 12 and D 12 at fixed Tffi. The last row 
gives the curvature coefficients obtained from Eq. (6] In 
the fits none of the corrections discussed in the text have 
been taken into account. 


fluctuations have been measured. As discussed earlier 
we consider E^^ to be a good approximation for the ra¬ 
tio E^-®. The difference between and can be 
corrected for by using the HRG model motivated correc¬ 
tion, R ^2 = This is appropriate for our NLO 

approximation and simply amounts to a rescaling of the 
abscissa in Fig. [21 While the intercept 012 in quadratic 
fits is not influenced by such a rescaling, C 12 increases by 
a factor (i?f 2 ^)~^i he. by about 20 % in the relevant tem¬ 
perature range. From Eq. |5| it is obvious that this will 
decrease the estimate for Analyzing the uncorrected 
STAR data thus will put an upper bound on /t^. 

Results for C 12 from quadratic fits to the two STAR 
data sets are given in Table HI We notice first that C 12 
extracted from the published STAR data 14, is about 
a factor 2 larger than the lattice QGD result for 05*2 (T/^o)- 
Since D 12 is positive, this corresponds to negative values 
for AC 2 as observed also in a HRG model analysis [ 2 ^, i.e. 
the STAR data on proton fluctuations taken in the range 
0.4GeV < pt < O.SGeV 141 are only compatible with 
a negative curvature coefficient k^- However, the still 
preliminary STAR data taken in the larger pt interval 
[l5j| . have a much smaller slope which is consistent with 
Ci 2 {Tf^o) within statistical errors. This gives an upper 
bound on the curvature of the freeze-out line. 


4 < 0.011 . 


( 8 ) 


Taking into account the HRG motivated correction for re¬ 
placing i?f 2 by Ri 2 reduces 4 makes the estimate for 
the upper bound compatible with zero, 4 = —0.012(15) 
for the STAR2.0 data set, i.e. the existing data for E^^ 
favor a small or even vanishing curvature of the freeze-out 
line at large ■s/s'^ as it is the case for the phenomeno¬ 
logical parametrization given in Ref. Q . 

f 

Let us finally compare the result obtained for a «2 with 
the curvature coefficient > 0 of the QGD transition 
line. 


Tc{pb) — Tcfi (1 — 4+ 0{4b)) ^ 


(9) 


















5 


where denotes the transition temperature at = 0 
Si- The above bound on is consistent with deter¬ 
minations of K 2 based on expansions of Tc{iib) around 
fJ.B = (> where the curvature coefficient is determined 
as a Taylor expansion coefficient evaluated at hb = 0. 
This gave ~ 0.007 Hi , which is about a factor 
two smaller than recent results for 3-261 that are 
based on lattice QCD calculations performed with large 
non-zero imaginary chemical potentials, corresponding 
to ^b/Tc^ — (1 — 3). These calculations yield values 
Kf ~ (0.015-0.02). 


CONCLUSIONS 

We provided a framework that allows to determine the 
curvature of the freeze-out line through a direct compar¬ 
ison between experimental data for mean and variance of 
net-electric charge and net-proton number fluctuations 
with lattice QCD calculations of cumulant ratios. We 
found the curvature of the freeze-out line to be small. It is 
consistent with the curvature of the QCD crossover line. 
At least for beam energies > 27 GeV our study 

suggests that freeze-out happens close to the crossover 
transition line. 

We have addressed some difficulties that arise when 
comparing lattice QCD calculations of conserved charge 
fluctuations with experimental data on cumulants of net- 
proton and net-charge fluctuations, although at present it 
is difficulty to provide a complete quantitative approach 
for this. We have pointed out some obvious differences 
between net-proton and net-baryon number fluctuations 
that are present already in equilibrium thermodynamics. 
We also addressed the question on the influence of trans¬ 
verse momentum cuts on the experimental data. Clearly 
control over these effects needs to be improved in fu¬ 
ture work. It also is known that cumulants of conserved 
charge fluctuations are sensitive to the width of the ra¬ 
pidity window covered by the experiments [ttI . [3 | . This 
dependence is more significant for net-electric charge fluc¬ 
tuations than for net-proton number fluctuations. In¬ 
creasing the rapidity window will decrease aq and thus 
will lead to an increase of E^. This will lead to larger val¬ 
ues for Tf^Q. However, at present these effects are difficult 
to quantify. 

Finally we note that throughout our analysis we as¬ 
sumed that the net-charge to net-baryon number ratio, 
r, is determined by the corresponding ratio present in the 
incident beams, r ~ 0.4. This ratio clearly will fluctu¬ 
ate in a fixed acceptance range covered in an experiment 
and one may argue that this ratio also shifts towards 
the isospin symmetric limit r = 0.5 for high beam ener¬ 
gies. Clearly, at present beam energies r < 0.5, as the 
net-charge expectation value would vanish in the isospin 
symmetric limit. In a more refined analysis one may, 
however, take into account effects arising from changes in 


the value of r which may also be beam energy dependent. 
In general, increasing r results in a decrease of £md 

f 

an increase of K 2 ■ Compared to other uncertainties, how¬ 
ever, this effect is small. For the preliminary STAR data 
set (STAR2.0) we find, for instance, that using r = 0.45 
decreases the by 3 MeV and shifts the curvature 
coefficient to a slightly positive value, = 0.004. 
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